Juvenile Pacific cod were tagged, acclimated to laboratory conditions, then gradually acclimated to four experimental temperatures (0, 5, 9, 16). After 6 weeks in treatments, fish were sacrificed, liver and fin tissues were collected for RNASeq and lcWGS, respectively, and measurements taken of fish length, wet weight, and liver weight.
On 11/21/22 160 tagged fish entered experimental tanks for acclimation, 40 per treatment (n=10 fish / tank, 4 tanks/treatment). On 12/28/22 temperatures began to slowly increase until all target experimental temperatures (0C, 5C, 9C, 16C) were reached on 1/8/23. Of the 40 fish per treatment, all survived the 0C, 5C, and 9C treatments. Four fish died in the 16degC treatment on 1/22/23 (tank 1), 1/24/23 (tank 2), 1/29/23 (tank 4), and 2/1/23 (tank 2). All survivors were euthanized/sampled on 2/8 & 2/9. Standard length and wet weight were 1) collected before tank acclimation, 2) before experimental temperature, and 3) at treatment termination so that growth and condition could be assessed prior to and during experimental treatments. Liver wet weight and tissues for sequencing were collected at treatment termination.
In this notebook I look at effects of temperature on distributions of growth rate (length, weight) body condition index (wet weight / length), and liver condition (i.e. hepato-somatic index, liver wet weight / whole body wet weight).
list.of.packages <- c("tidyverse", "readxl", "janitor", "purrr", "ggpubr", "googlesheets4", "plotly", "lubridate")
# Load all required libraries
all.packages <- c(list.of.packages)
lapply(all.packages, FUN = function(X) {
do.call("require", list(X))
})
`%!in%` = Negate(`%in%`)
load(file = "../rnaseq/aligned/counts.ts")
phenotypes <- read_excel("../data/Pcod Temp Growth experiment 2022-23 DATA.xlsx", sheet = "AllData") %>% clean_names() %>%
mutate_at(c("tank", "temperature", "microchip_id", "dissection_date", "genetic_sampling_count"), factor) %>%
dplyr::rename("sl_final"="sl_mm", "wwt_final"="whole_body_ww_g") %>%
mutate(growth.sl.accl=sl_12272022-sl_11212022,
growth.sl.trt=round(sl_final-sl_12272022, digits = 2),
growth.wwt.trt=wwt_final-wwt_12272022,
growth.wwt.accl=round(wwt_12272022-wwt_11212022, digits = 2),
ci=wwt_final/sl_final,
hsi=total_liver_ww_mg/(wwt_final*1000),
mort=as.factor(case_when(is.na(mort_date) ~ "No", TRUE ~ "Yes")),
rna=case_when(genetic_sampling_count %in% gsub("sample_", "", rownames(counts.ts)) ~ "Yes", TRUE ~ "No"),
measure_date=case_when(is.na(mort_date) ~ ymd(dissection_date), TRUE ~ mort_date)) %>%
mutate(days_growth = as.numeric(difftime(ymd(measure_date), ymd("2022-12-27"), units = "days")))
phenotypes %>% head()
## # A tibble: 6 × 32
## microchip_id sl_11212022 wwt_11212022 tank temperature sl_12272022
## <fct> <dbl> <dbl> <fct> <fct> <dbl>
## 1 9397 86 7.38 12 0 93
## 2 9467 92 7.98 12 0 97
## 3 5059 102 11.6 12 0 113
## 4 9461 82 6.15 12 0 88
## 5 9560 95 9.59 12 0 101
## 6 7386 95 9.23 12 0 108
## # ℹ 26 more variables: wwt_12272022 <dbl>, mort_date <dttm>,
## # dissection_date <fct>, sl_final <dbl>, wwt_final <dbl>,
## # total_liver_ww_mg <dbl>, liverfor_lipids_ww_mg <dbl>,
## # muscle_w_wfor_lipids_mg <dbl>, genetic_sampling_count <fct>,
## # dissection_comments <chr>, molecular_notes <chr>, fin <chr>, liver <chr>,
## # blood <chr>, spleen <chr>, gill <chr>, growth.sl.accl <dbl>,
## # growth.sl.trt <dbl>, growth.wwt.trt <dbl>, growth.wwt.accl <dbl>, …
# Wet weight over time
phenotypes %>%
group_by(temperature) %>%
dplyr::summarise(sl_mean.1=mean(sl_11212022),sl_sd.1=sd(sl_11212022),
sl_mean.2=mean(sl_12272022),sl_sd.2=sd(sl_12272022),
sl_mean.3=mean(sl_final),sl_sd.3=sd(sl_final),
wwt_mean.1=mean(wwt_11212022),wwt_sd.1=sd(wwt_11212022),
wwt_mean.2=mean(wwt_12272022),wwt_sd.2=sd(wwt_12272022),
wwt_mean.3=mean(wwt_final),wwt_sd.3=sd(wwt_final)) %>%
pivot_longer(cols = -temperature) %>%
separate(name, sep = "\\.", into = c("metric", "time")) %>%
mutate_at(c("metric"), factor) %>% mutate(time=as.numeric(time)) %>%
filter(metric=="wwt_mean") %>%
ggplot() + geom_line(aes(x=time, y=value, color=temperature), cex=1.5) + theme_minimal() +
scale_color_manual(name="Temperature",
values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
ggtitle("Wet weight over time") + xlab(NULL) + ylab("Wet weight") +
scale_x_continuous(breaks = c(1, 2, 3), labels = c("Pre-Acclimation", "Pre-Treatment", "Treatment\ntermination"))
### Standard length by temperature over time - mean and SD
phenotypes %>%
group_by(temperature) %>%
dplyr::summarise(sl_mean.1=mean(sl_11212022),sl_sd.1=sd(sl_11212022),
sl_mean.2=mean(sl_12272022),sl_sd.2=sd(sl_12272022),
sl_mean.3=mean(sl_final),sl_sd.3=sd(sl_final),
wwt_mean.1=mean(wwt_11212022),wwt_sd.1=sd(wwt_11212022),
wwt_mean.2=mean(wwt_12272022),wwt_sd.2=sd(wwt_12272022),
wwt_mean.3=mean(wwt_final),wwt_sd.3=sd(wwt_final)) %>%
pivot_longer(cols = -temperature) %>%
separate(name, sep = "\\.", into = c("metric", "time")) %>%
mutate_at(c("metric"), factor) %>% mutate(time=as.numeric(time)) %>%
filter(metric=="sl_mean") %>%
ggplot() + geom_line(aes(x=time, y=value, color=temperature), cex=1.5) + theme_minimal() +
scale_color_manual(name="Temperature",
values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
ggtitle("Standard length over time") + xlab(NULL) + ylab("Standard length") +
scale_x_continuous(breaks = c(1, 2, 3), labels = c("Pre-Acclimation", "Pre-Treatment", "Treatment\ntermination"))
phenotypes %>%
ggplot(aes(y=growth.sl.accl/days_growth, x=temperature, fill=temperature, shape=mort, alpha=mort)) +
geom_violin(trim = F) + geom_point(position = position_jitterdodge(jitter.width = 1)) + theme_minimal() +
ggtitle("Growth rate while acclimating, standard length (mm/day)") +
scale_alpha_manual(values=c(0.75,0.5)) +
scale_x_discrete(drop=T) + #Do drop empty factors for temperature contrasts
ylab("mm / day") +
scale_fill_manual(name="Temperature", values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
xlab(NULL) + ylab(NULL)
# phenotypes %>% filter(mort=="No") %>%
# ggplot(aes(y=growth.sl.accl/days_growth, x=temperature, fill=temperature)) +
# geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
# ggtitle("Growth rate while acclimating, standard length (mm/day)") +
# scale_fill_manual(name="Temperature",
# values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
# xlab(NULL) + ylab(NULL)
phenotypes %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(growth.sl.accl), sd=sd(growth.sl.accl))
## # A tibble: 4 × 3
## temperature mean sd
## <fct> <dbl> <dbl>
## 1 0 8.7 2.82
## 2 5 8.43 3.34
## 3 9 8 2.42
## 4 16 9.78 2.75
summary(aov(growth.sl.accl ~ temperature, phenotypes %>% filter(mort=="No")))
## Df Sum Sq Mean Sq F value Pr(>F)
## temperature 3 64.4 21.481 2.637 0.0518 .
## Residuals 152 1238.4 8.147
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov(growth.sl.accl ~ temperature, phenotypes %>% filter(mort=="No")))
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = growth.sl.accl ~ temperature, data = phenotypes %>% filter(mort == "No"))
##
## $temperature
## diff lwr upr p adj
## 5-0 -0.275000 -1.93297778 1.3829778 0.9730922
## 9-0 -0.700000 -2.35797778 0.9579778 0.6920681
## 16-0 1.077778 -0.62563246 2.7811880 0.3574449
## 9-5 -0.425000 -2.08297778 1.2329778 0.9097552
## 16-5 1.352778 -0.35063246 3.0561880 0.1700574
## 16-9 1.777778 0.07436754 3.4811880 0.0371523
phenotypes %>%
ggplot(aes(y=growth.sl.trt/days_growth, x=temperature, fill=temperature, shape=mort, alpha=mort)) +
geom_violin() + geom_point(position = position_jitterdodge(jitter.width = 1)) + theme_minimal() +
ggtitle("Growth rate in treatments, standard length (mm/day)") +
scale_alpha_manual(values=c(0.75,0.5)) +
ylab("mm / day") +
scale_fill_manual(name="Temperature", values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
xlab(NULL) + ylab(NULL)
# phenotypes %>% filter(mort=="No") %>%
# ggplot(aes(y=growth.sl.trt, x=temperature, fill=temperature)) +
# geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
# ggtitle("Growth in treatments, standard length") +
# scale_fill_manual(name="Temperature",
# values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c"))
summary(aov(growth.sl.trt ~ temperature, phenotypes %>% filter(mort=="No")))
## Df Sum Sq Mean Sq F value Pr(>F)
## temperature 3 2456 818.8 87.06 <2e-16 ***
## Residuals 152 1430 9.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov(growth.sl.trt ~ temperature, phenotypes %>% filter(mort=="No")))
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = growth.sl.trt ~ temperature, data = phenotypes %>% filter(mort == "No"))
##
## $temperature
## diff lwr upr p adj
## 5-0 4.875000 3.093614 6.6563863 0.0000000
## 9-0 10.275000 8.493614 12.0563863 0.0000000
## 16-0 8.569444 6.739244 10.3996449 0.0000000
## 9-5 5.400000 3.618614 7.1813863 0.0000000
## 16-5 3.694444 1.864244 5.5246449 0.0000031
## 16-9 -1.705556 -3.535756 0.1246449 0.0774644
phenotypes %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(growth.sl.trt), sd=sd(growth.sl.trt))
## # A tibble: 4 × 3
## temperature mean sd
## <fct> <dbl> <dbl>
## 1 0 4.12 1.92
## 2 5 9 2.53
## 3 9 14.4 3.51
## 4 16 12.7 3.98
phenotypes %>%
ggplot(aes(y=growth.wwt.accl/days_growth, x=temperature, fill=temperature, shape=mort, alpha=mort)) +
geom_violin() + geom_point(position = position_jitterdodge(jitter.width = 1)) + theme_minimal() +
ggtitle("Growth rate while acclimating, wet weight (g/day)") +
scale_alpha_manual(values=c(0.75,0.5)) +
ylab("g / day") +
scale_fill_manual(name="Temperature", values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
xlab(NULL)
# phenotypes %>% filter(mort=="No") %>%
# ggplot(aes(y=growth.wwt.accl, x=temperature, fill=temperature)) +
# geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
# ggtitle("Growth while acclimating, wet weight") +
# scale_fill_manual(name="Temperature",
# values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c"))
summary(aov(growth.wwt.accl ~ temperature, phenotypes %>% filter(mort=="No")))
## Df Sum Sq Mean Sq F value Pr(>F)
## temperature 3 3.64 1.2122 1.85 0.14
## Residuals 152 99.59 0.6552
TukeyHSD(aov(growth.wwt.accl ~ temperature, phenotypes %>% filter(mort=="No")))
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = growth.wwt.accl ~ temperature, data = phenotypes %>% filter(mort == "No"))
##
## $temperature
## diff lwr upr p adj
## 5-0 -0.23050000 -0.7006701 0.2396701 0.5811808
## 9-0 -0.32350000 -0.7936701 0.1466701 0.2834053
## 16-0 0.04158333 -0.4414705 0.5246372 0.9960401
## 9-5 -0.09300000 -0.5631701 0.3771701 0.9556859
## 16-5 0.27208333 -0.2109705 0.7551372 0.4623147
## 16-9 0.36508333 -0.1179705 0.8481372 0.2065730
phenotypes %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(growth.wwt.accl), sd=sd(growth.wwt.accl))
## # A tibble: 4 × 3
## temperature mean sd
## <fct> <dbl> <dbl>
## 1 0 2.45 0.994
## 2 5 2.22 0.682
## 3 9 2.13 0.811
## 4 16 2.50 0.703
phenotypes %>%
ggplot(aes(y=growth.wwt.trt/days_growth, x=temperature, fill=temperature, shape=mort, alpha=mort)) +
geom_violin() + geom_point(position = position_jitterdodge(jitter.width = 1)) + theme_minimal() +
ggtitle("Growth rate in treatments, wet weight (g/day)") +
scale_alpha_manual(values=c(0.75,0.5)) +
ylab("g / day") +
scale_fill_manual(name="Temperature", values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c"))
# phenotypes %>% filter(mort=="No") %>%
# ggplot(aes(y=growth.wwt.trt, x=temperature, fill=temperature)) +
# geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
# ggtitle("Growth in treatments, wet weight") +
# scale_fill_manual(name="Temperature",
# values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c"))
summary(aov(growth.wwt.trt ~ temperature, phenotypes %>% filter(mort=="No")))
## Df Sum Sq Mean Sq F value Pr(>F)
## temperature 3 283.3 94.44 92.34 <2e-16 ***
## Residuals 152 155.4 1.02
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov(growth.wwt.trt ~ temperature, phenotypes %>% filter(mort=="No")))
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = growth.wwt.trt ~ temperature, data = phenotypes %>% filter(mort == "No"))
##
## $temperature
## diff lwr upr p adj
## 5-0 1.4110000 0.8235964 1.9984036 0.0000000
## 9-0 3.3832500 2.7958464 3.9706536 0.0000000
## 16-0 2.9721389 2.3686391 3.5756387 0.0000000
## 9-5 1.9722500 1.3848464 2.5596536 0.0000000
## 16-5 1.5611389 0.9576391 2.1646387 0.0000000
## 16-9 -0.4111111 -1.0146109 0.1923887 0.2920273
phenotypes %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(growth.wwt.trt), sd=sd(growth.wwt.trt))
## # A tibble: 4 × 3
## temperature mean sd
## <fct> <dbl> <dbl>
## 1 0 0.0518 0.526
## 2 5 1.46 0.770
## 3 9 3.44 1.25
## 4 16 3.02 1.32
phenotypes %>%
ggplot(aes(y=ci, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) +
geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
ggtitle("Body condition index\n(wet weight / standard length") +
scale_alpha_manual(values=c(0.75,0.5)) +
scale_fill_manual(name="Temperature",
values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
xlab(NULL) + ylab("Wet weight / standard length")
# phenotypes %>% filter(mort=="No") %>%
# ggplot(aes(y=ci, x=temperature, fill=temperature)) +
# geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
# ggtitle("Appr. body condition index\n(wet weight / standard length") +
# scale_fill_manual(name="Temperature",
# values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c"))
summary(aov(ci ~ temperature, phenotypes %>% filter(mort=="No")))
## Df Sum Sq Mean Sq F value Pr(>F)
## temperature 3 0.00965 0.003216 10.39 2.91e-06 ***
## Residuals 152 0.04702 0.000309
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov(ci ~ temperature, phenotypes %>% filter(mort=="No")))
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = ci ~ temperature, data = phenotypes %>% filter(mort == "No"))
##
## $temperature
## diff lwr upr p adj
## 5-0 0.002714042 -0.007502273 0.01293036 0.9007512
## 9-0 0.014829986 0.004613671 0.02504630 0.0013176
## 16-0 0.018819590 0.008323324 0.02931586 0.0000408
## 9-5 0.012115944 0.001899629 0.02233226 0.0129654
## 16-5 0.016105547 0.005609282 0.02660181 0.0005974
## 16-9 0.003989604 -0.006506662 0.01448587 0.7568554
phenotypes %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(ci), sd=sd(ci))
## # A tibble: 4 × 3
## temperature mean sd
## <fct> <dbl> <dbl>
## 1 0 0.104 0.0168
## 2 5 0.107 0.0140
## 3 9 0.119 0.0214
## 4 16 0.123 0.0174
phenotypes %>%
ggplot(aes(y=hsi, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort)) +
geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
ggtitle("Hepato-somatic index\n(liver wet weight / total wet weight)") +
scale_alpha_manual(values=c(0.75,0.5)) +
scale_fill_manual(name="Temperature",
values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
xlab(NULL) + ylab("Liver wet weight / Total wet weight")
# phenotypes %>% filter(mort=="No") %>%
# ggplot(aes(y=hsi, x=temperature, fill=temperature)) +
# geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
# ggtitle("Appr. hepato-somatic index\n(liver wet weight / total wet weight)") +
# scale_fill_manual(name="Temperature",
# values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c"))
summary(aov(hsi ~ temperature, phenotypes %>% filter(mort=="No")))
## Df Sum Sq Mean Sq F value Pr(>F)
## temperature 3 4.945e-09 1.648e-09 31.45 7.14e-16 ***
## Residuals 152 7.966e-09 5.240e-11
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov(hsi ~ temperature, phenotypes %>% filter(mort=="No")))
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = hsi ~ temperature, data = phenotypes %>% filter(mort == "No"))
##
## $temperature
## diff lwr upr p adj
## 5-0 -2.335143e-06 -6.540247e-06 1.869961e-06 0.4749567
## 9-0 -8.358065e-06 -1.256317e-05 -4.152961e-06 0.0000045
## 16-0 -1.473728e-05 -1.905762e-05 -1.041695e-05 0.0000000
## 9-5 -6.022922e-06 -1.022803e-05 -1.817818e-06 0.0015764
## 16-5 -1.240214e-05 -1.672247e-05 -8.081807e-06 0.0000000
## 16-9 -6.379219e-06 -1.069955e-05 -2.058885e-06 0.0010420
phenotypes %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(hsi), sd=sd(hsi))
## # A tibble: 4 × 3
## temperature mean sd
## <fct> <dbl> <dbl>
## 1 0 0.0000361 0.00000674
## 2 5 0.0000337 0.00000877
## 3 9 0.0000277 0.00000704
## 4 16 0.0000213 0.00000600
Within treatments, I’d like to use WGCNA to identify co-expressed genes that are linearly related to a performance indicator. Here, I explore using HSI*CI as the performance indicator. In 16C, HSI was lower than control (9C) but CI was higher than control, so fish with high CI & HSI could be “high performers”, and expression patterns could provide molecular indicators of performance.
# Condition Index * Hepato-somatic Index
phenotypes %>%
ggplot(aes(y=ci*hsi, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) +
geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
ggtitle("Possible index of performance (?)\n(Condition index * Hepatosomatic Index)") +
scale_alpha_manual(values=c(0.75,0.5)) +
scale_fill_manual(name="Temperature",
values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
xlab(NULL) + ylab("CI * HSI")
summary(aov(hsi*ci ~ temperature, phenotypes %>% filter(mort=="No")))
## Df Sum Sq Mean Sq F value Pr(>F)
## temperature 3 2.921e-11 9.738e-12 9.545 8.18e-06 ***
## Residuals 152 1.551e-10 1.020e-12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov(hsi*ci ~ temperature, phenotypes %>% filter(mort=="No")))
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = hsi * ci ~ temperature, data = phenotypes %>% filter(mort == "No"))
##
## $temperature
## diff lwr upr p adj
## 5-0 -1.648471e-07 -7.515427e-07 4.218485e-07 0.8849743
## 9-0 -4.366931e-07 -1.023389e-06 1.500025e-07 0.2184555
## 16-0 -1.156764e-06 -1.759537e-06 -5.539917e-07 0.0000099
## 9-5 -2.718460e-07 -8.585416e-07 3.148496e-07 0.6254920
## 16-5 -9.919170e-07 -1.594689e-06 -3.891446e-07 0.0001958
## 16-9 -7.200711e-07 -1.322843e-06 -1.172986e-07 0.0121118
phenotypes %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(hsi*ci), sd=sd(hsi*ci))
## # A tibble: 4 × 3
## temperature mean sd
## <fct> <dbl> <dbl>
## 1 0 0.00000378 0.00000101
## 2 5 0.00000362 0.00000106
## 3 9 0.00000335 0.00000110
## 4 16 0.00000263 0.000000828
#ggplotly(
phenotypes %>% filter(mort=="No") %>%
ggplot(aes(y=hsi, x=ci, color=temperature)) +
geom_point() + theme_minimal() +
ggtitle("HSI / CI") +
scale_color_manual(name="Temperature",
values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c"))+
geom_smooth(method = "lm")#)
ggplotly(phenotypes %>% filter(mort=="No") %>%
ggplot(aes(y=hsi*ci, x=interaction(rna, temperature), fill=temperature, label=genetic_sampling_count, shape=rna)) +
geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
ggtitle("Performance Index (CI*HSI) by temperature, by RNASeq data or none") +
scale_fill_manual(name="Temperature",
values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")))
# ggplotly(phenotypes %>% filter(rna=="Yes") %>%
# ggplot(aes(y=hsi*ci, x=temperature, fill=temperature, label=genetic_sampling_count)) +
# geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
# ggtitle("Appr. hepato-somatic index\n(liver wet weight / total wet weight)") +
# scale_alpha_manual(values=c(0.75,0.5)) + ggtitle("Performance index by temperature, RNASeq samples only") +
# scale_fill_manual(name="Temperature",
# values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")))
#
# ggplotly(phenotypes %>% filter(rna=="Yes") %>%
# left_join(haplos.allzp3, by = c("genetic_sampling_count"="sample")) %>%
# ggplot(aes(y=hsi*ci, x=temperature, fill=temperature, label=genetic_sampling_count, shape=haplo.exons.lcwgs)) +
# geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
# ggtitle("Appr. hepato-somatic index\n(liver wet weight / total wet weight)") +
# scale_fill_manual(name="Temperature",
# values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")))
lipids <- read_excel("../data/Lipid class liver data_091324.xlsx") %>% clean_names() %>%
mutate(sample_id=as.factor(sample_id))
## New names:
## • `TAG` -> `TAG...20`
## • `FFA` -> `FFA...21`
## • `ST` -> `ST...22`
## • `Polar Lipids` -> `Polar Lipids...23`
## • `Total` -> `Total...24`
## • `TAG` -> `TAG...28`
## • `FFA` -> `FFA...29`
## • `ST` -> `ST...30`
## • `Polar Lipids` -> `Polar Lipids...31`
## • `Total` -> `Total...32`
phenotypes2 <- phenotypes %>%
left_join(lipids %>% dplyr::select(sample_id, genetics_sample_number, tag_20:total_32),
by=c("microchip_id"="sample_id")) %>%
mutate(pi=1000*growth.sl.trt*ci*hsi*total_32)
phenotypes2 %>%
filter(!is.na(total_32)) %>%
ggplot(aes(y=tag_28, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) +
geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
ggtitle("Triglyceride content (ug/mg)") +
scale_alpha_manual(values=c(0.75,0.5)) +
scale_fill_manual(name="Temperature",
values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
xlab(NULL) + ylab("TAG Content")
## Warning: Groups with fewer than two data points have been dropped.
## Warning: The following aesthetics were dropped during statistical transformation: label
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
phenotypes2 %>%
filter(!is.na(total_32)) %>%
ggplot(aes(y=ffa_29, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) +
geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
ggtitle("Free fatty acid content (ug/mg)") +
scale_alpha_manual(values=c(0.75,0.5)) +
scale_fill_manual(name="Temperature",
values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
xlab(NULL) + ylab("FFA Content")
## Warning: Groups with fewer than two data points have been dropped.
## The following aesthetics were dropped during statistical transformation: label
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
phenotypes2 %>%
filter(!is.na(total_32)) %>%
ggplot(aes(y=st_30, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) +
geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
ggtitle("Sterol content (ug/mg)") +
scale_alpha_manual(values=c(0.75,0.5)) +
scale_fill_manual(name="Temperature",
values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
xlab(NULL) + ylab("sterol Content")
## Warning: Groups with fewer than two data points have been dropped.
## The following aesthetics were dropped during statistical transformation: label
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
phenotypes2 %>%
filter(!is.na(total_32)) %>%
ggplot(aes(y=polar_lipids_31, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) +
geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
ggtitle("Polar lipid content (ug/mg)") +
scale_alpha_manual(values=c(0.75,0.5)) +
scale_fill_manual(name="Temperature",
values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
xlab(NULL) + ylab("Polar lipid Content")
## Warning: Groups with fewer than two data points have been dropped.
## The following aesthetics were dropped during statistical transformation: label
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
phenotypes2 %>%
filter(!is.na(total_32)) %>%
ggplot(aes(y=total_32, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) +
geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
ggtitle("Total lipid content (ug/mg)") +
scale_alpha_manual(values=c(0.75,0.5)) +
scale_fill_manual(name="Temperature",
values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
xlab(NULL) + ylab("Total lipid Content")
## Warning: Groups with fewer than two data points have been dropped.
## The following aesthetics were dropped during statistical transformation: label
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
# Condition Index * Hepato-somatic Index * Total Lipid
#ggplotly(
phenotypes2 %>%
# First I scale each metric so they are between about 1-20
mutate(hsi=hsi*1e5,
ci=ci*100,
total_32=total_32/100) %>%
#then I calculate a performance index from all metrics (x100)
mutate(pi=(growth.sl.trt*hsi*ci*total_32)/100) %>%
ggplot(aes(y=pi, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) +
geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
ggtitle("Possible composite index of performance\n(Growth * Condition index * Hepatosomatic Index * Total Lipid)") +
scale_alpha_manual(values=c(0.75,0.5)) +
scale_fill_manual(name="Temperature",
values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
xlab(NULL) + ylab("Growth + CI * HSI * Total Lipid * 1,000")
#)